III. Soving the diffusion problem with implicit Euler integration with multigrid & finite differences


In [1]:
using PyPlot


INFO: Loading help data...

Problem description

For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.

\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}

We therefore want to solve the following PDE:

$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$

For this third solution, we again use backward Euler integration in time and a second-order finite differences method in space. This leads to the following method to calculate the next time step:

$$ u^{n+1} = u^n + D \cdot \Delta t \cdot A u^{n+1} $$

Here, A is the finite differences matrix for the 2D Laplace operator. However, as $u^{n+1}$ is not known in advance, we have to solve this equation for it. The scheme becomes the following, which has to be solved for $u^{n+1}$:

$$ \left( I - D \cdot \Delta t \cdot A \right) u^{n+1} = u^{n} $$

In this case, however, we don’t solve for $u^{n+1}$ directly (e.g. with $LU$-decomposition). Instead, we use a multigrid approach to solve the system iteratively. As the spatial approximation with finite differences has a simple, regular form, we can easily create the finite differences matrix for other, coarser grids as well. The matrices have the exact same form as for the finest grid, except that they have fewer entries. However, we impose that the number of intervals is a power of two to ensure that they can easily be mapped to a coarser grid.

Types for spatial and temporal domain


In [2]:
immutable Grid2D
    
    x_start::Real
    x_end::Real
    y_start::Real
    y_end::Real
    Nx::Int
    Ny::Int
    dx::Real
    dy::Real
    levels::Int
    
    function Grid2D(x_start, x_end, y_start, y_end, Nx, Ny, levels)
        
        if Nx/2^(levels-1) != int(Nx/2^(levels-1))
            error("ERROR: The number of intervals in x-direction is not divisible by 2^i")
        end
        
        if Ny/2^(levels-1) != int(Ny/2^(levels-1))
            error("ERROR: The number of intervals in y-direction is not divisible by 2^i")
        end
        
        new(x_start, x_end, y_start, y_end, Nx, Ny, (x_end-x_start)/Nx, (y_end-y_start)/Ny, levels)
    end
end

In [3]:
immutable TimeSteps
    
    t_start::Real
    t_end::Real
    Nt::Int
    dt::Real
    
    TimeSteps(t_start, t_end, Nt) = new(t_start, t_end, Nt, (t_end-t_start)/Nt)
    
end

Generator for exact solution

As a very simple scenario, we use the following initial conditions:

$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$

This leads to the following analytical solution:

$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$

In [4]:
function gen_exact(grid, C0, D)
    x_size = grid.x_end - grid.x_start
    y_size = grid.y_end - grid.y_start
    exact(x,y,t) = C0 * sin((x-grid.x_start)/x_size * pi) * sin((y-grid.y_start)/y_size * pi) *
                   exp( - (1/x_size^2 + 1/y_size^2) * pi^2 * D * t)
end


Out[4]:
gen_exact (generic function with 1 method)

Matrix for finite differences

We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).

Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:

\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}

The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.

\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}

The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.

$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$

As we have to create several different matrices for the multigrid matrix, we define a function that returns those matrices depending on the number and size of intervals.


In [5]:
function fdmatrix(grid::Grid2D, matrix_function::Function = A -> A)
    
    # number of unknowns in x and y direction
    ux = (grid.Nx-1)
    uy = (grid.Ny-1)
    
    # finite difference matrix for ∂²u/∂x²
    FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1))
    FDx = blkdiag({FDx_local for i in 1:uy}...)
    
    # finite difference matrix for ∂²u/∂y²
    FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux))
    
    # return combined matrix for Δu
    matrix_function(1/grid.dx^2 * FDx + 1/grid.dy^2 * FDy)
end


Out[5]:
fdmatrix (generic function with 2 methods)

Multigrid for spatial derivatives


In [6]:
# helper function to get coarser levels of a grid
function level(grid::Grid2D, lvl)
    if lvl > grid.levels
        error("ERROR: Grid doesn't support level $(level).")
    end
    Grid2D(grid.x_start, grid.x_end, grid.y_start, grid.y_end, int(grid.Nx/2^(lvl-1)), int(grid.Ny/2^(lvl-1)), grid.levels-(lvl-1))
end


Out[6]:
level (generic function with 1 method)

In [7]:
function weighted_jacobi(A, x, b)
    weight = 2/3
    (speye(A) - weight*A) * x + weight*b
end


Out[7]:
weighted_jacobi (generic function with 1 method)

In [8]:
function restrictionmatrix(grid::Grid2D)
    ux_old = grid.Nx - 1
    uy_old = grid.Ny - 1
    ux_new = div(grid.Nx,2) - 1
    uy_new = div(grid.Ny,2) - 1
    R = spzeros(ux_new*uy_new, ux_old*uy_old)
    for j = 1:uy_new
        for i = 1:ux_new
            i_new = ux_new*(j-1) + i
            R[i_new, ux_old*(2*j-1)     + 2*i    ] = 4
            R[i_new, ux_old*(2*j-1)     + 2*i - 1] = 2
            R[i_new, ux_old*(2*j-1)     + 2*i + 1] = 2
            R[i_new, ux_old*(2*j-1 - 1) + 2*i    ] = 2
            R[i_new, ux_old*(2*j-1 - 1) + 2*i - 1] = 1
            R[i_new, ux_old*(2*j-1 - 1) + 2*i + 1] = 1
            R[i_new, ux_old*(2*j-1 + 1) + 2*i    ] = 2
            R[i_new, ux_old*(2*j-1 + 1) + 2*i - 1] = 1
            R[i_new, ux_old*(2*j-1 + 1) + 2*i + 1] = 1
        end
    end
    R /= 16
end


Out[8]:
restrictionmatrix (generic function with 1 method)

In [9]:
function gen_mgsolver(matrix_generator, restriction_generator, smoother, grid::Grid2D, ν1::Int, ν2::Int, tol=1e-6, max_it = 100)
    
    # create main matrices for all levels
    A = Dict()
    for i = 1:grid.levels
        A[i] = matrix_generator(level(grid,i))
    end
    
    # create restriction matrices for all levels
    R = Dict()
    for i = 1:grid.levels-1
        R[i] = restriction_generator(level(grid,i))
    end
    
    function mgsolver(b::Array, x0::Array = b, level::Int = 1)
        
        if level == grid.levels
            return A[level] \ b
        end
        
        x = x0
        
        for it = 1:max_it
                
            # smoothen problem on current level
            for i = 1:ν1
                x = smoother(A[level], x, b)
            end
            
            # restrict residual to next level and solve recursively
            r = R[level] * (b - A[level] * x)
            e = mgsolver(r, zeros(r), level+1)
        
            # prolongate error to current level and add to solution
            x += 4*R[level].' * e
        
            # post-smoothening on current level
            for i = 1:ν2
                x = smoother(A[level], x, b)
            end
            
            if level > 1
                return x
            end
            
            if norm(A[level] * x - b) < tol
                #println("Converged with $(it) iterations")
                return x
            end
        end
    
        return x
    end
    
    return mgsolver
end


Out[9]:
gen_mgsolver (generic function with 3 methods)

Do integration with backward Euler method


In [10]:
function integrate(x::Array, integration_step::Function, time::TimeSteps)
    for i in 1:time.Nt
        x = integration_step(x)
    end
    x
end


Out[10]:
integrate (generic function with 1 method)

Define and solve the problem


In [11]:
g = Grid2D(1,3,2,5,2^6,2^6,3)


Out[11]:
Grid2D(1,3,2,5,64,64,0.03125,0.046875,3)

In [12]:
t = TimeSteps(0, 20, 100)


Out[12]:
TimeSteps(0,20,100,0.2)

In [13]:
C0 = 5
D  = 0.001
exact = gen_exact(g, C0, D)


Out[13]:
exact (generic function with 1 method)

Exact solution


In [14]:
C_exact = Float64[exact(x,y,t.t_end) for x=g.x_start:g.dx:g.x_end, y=g.y_start:g.dy:g.y_end]
pcolormesh(C_exact)
colorbar()
clim((0,5))


Approximate solution


In [15]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,t.t_start) for x=g.x_start:g.dx:g.x_end, y=g.y_start:g.dy:g.y_end];

In [16]:
integration_step = gen_mgsolver(g -> speye((g.Nx-1)*(g.Ny-1)) - t.dt * D * fdmatrix(g),
                                restrictionmatrix, weighted_jacobi, g, 0, 1)


Out[16]:
mgsolver (generic function with 3 methods)

In [17]:
x = vcat(C[2:end-1,2:end-1]...)
x = integrate(x, integration_step, t)
C[2:end-1,2:end-1] = reshape(x,g.Nx-1,g.Ny-1)
pcolormesh(C)
colorbar()
clim((0,C0))


Compare to exact solution


In [18]:
vecnorm(C - C_exact)


Out[18]:
0.005914204557685913